 function [X,W,WInv,WChol]=MNIWDraw(muMat,PChol,SInv,v)
%  
% Generates a draw of (X,W)~MNIW(muMat,P,S,v)
% such that 
% X|W ~ MN( muMat, W   kron P ) 
%   W ~ IW( v    , S          ) 
% 
% *Notes* 
% 1. *Pchol=chol(P)*
% 2. *Sinv=inv(S)*
% 
%% Input 
% muMat: [p,q] matrix with Mean 
% Pprec: [p,p] matrix with *PChol=Chol(P)*, i.e. P=PChol'*PChol
% SInv:  [q,q] INVERSE matrix for IW, Sinv=Inv(S) 
% v:     (scalar) degrees of Freedom for IW 
%
%% Output 

% X  [p,q] draw from  X|W ~ MN( muMat, W   kron P ) 
% W  [q,q] draw from  IW( v , S )   S=inv(SInv)
% WInv  [q,q] inverse of W through SVD if Nargout > 2 
% WChol [q,q] chol    of W s.t. W=WChol'*WChol
[Nr,Nc]=size(muMat); 
% Nr=k 
% Nc=n 
%%
drMat=mvnrnd(zeros(Nc,1),SInv,v); 
W=eye(Nc)/(drMat'*drMat); 
%% 1. Obtain chol(W) and inv(W) using the SVD 
W=0.5*(W+W'); 
% This is more robust but probably slower than inv(W)
% for n small 
%% 1.a SVD 
% PP*DD*PPinv'=W  Notice the transpose 
% PPinv=inv(PP)' 
% PPinv'=inv(PP); 
[PP,DD,PPinv]=svd(W); 

%% 1.b Truncate small singular values 
tolZero=1e-12; 
firstZero=min(find(diag(DD) < tolZero));
if isempty(firstZero)==false 
    PP=PP(:,1:firstZero-1); 
    PPinv=PPinv(:,1:firstZero-1); 
    DD=DD(1:firstZero-1,1:firstZero-1); 
end 
%% 1.c Inverse and Cholesky 
WChol=sqrt(DD)*PPinv';
X=(PChol')*(randn(Nr,Nc))*WChol+muMat;
if nargout > 2
    WInv =PP*(DD\eye(Nc))*PPinv';
end